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Abstract. Recently, we proposed a self-propelled particle model with 
competing alignment interactions: nearby particles tend to align their 
velocities whereas they anti-align their direction of motion with par¬ 
ticles which are further away [R. Grofimann et al., Phys. Rev. Lett. 
113 , 258104 (2014)]. Here, we extend our previous numerical analysis 
of the high density regime considering low particle densities too. We 
report on the emergence of various macroscopic patterns such as vortex 
arrays, mesoscale turbulence as well as the formation of polar clusters, 
polar bands and nematically ordered states. Furthermore, we study 
analytically the instabilities leading to pattern formation in mean-field 
approximation. We argue that these instabilities are well described by 
a reduced set of hydrodynamic equations in the limit of high density. 


1 Introduction 


The study of active matter is one of the main topics in modern classical statistical 
mechanics [T|j5]. Active matter systems are composed of self-driven entities, so called 
active particles, which take up energy from their environment in order to move actively 
in a dissipative medium 2 . Due to the complex interplay of self-propulsion and 
dissipation, the dynamics of active matter systems is inherently out of thermodynamic 
equilibrium. Thus, complex spatio-temporal patterns may emerge in interacting active 
particle systems. There are numerous experimental realizations of active matter, for 
example driven filaments on the molecular scale [6] , colloidal systems 17,8], bacterial 
colonies (9-17 as well as macroscopic systems like flocking birds 18 and schooling 
fish fl9[ 


The experimental study of bacterial colonies has become increasingly important 
for the understanding of active matter. On the one hand, the quantitative analysis 
of fascinating phenomena like clustering 11 and vortex formation 20 provide an 


a e-mail: grossmann@physik.hu-berlin.de 
b e-mail: pawelr@princeton. edu 










2 


Will be inserted by the editor 


excellent opportunity to test theoretical predictions, and, on the other hand, exper¬ 
iments often reveal unforeseen novel phases of active matter. Recently, a number of 
experiments reported the emergence of dynamic vortex patterns in dense suspensions 
of B acil lus subtilis [12 -15], Paenibacillus dendritiformis 16 as well as Escherichia 
coli 17 . In particular, it was observed that bacteria locally align their body axis. 
However, dynamic vortices (defects) in the velocity field appeared at larger scales. 
Interestingly, vortices are separated by a characteristic length scale which is larger 
than the size of an individual bacterium but smaller than the system size. This novel 
phenomenon was termed mesoscale turbulence. 

Mesoscale turbulence was suggested to be described by a phenomenological hydro- 

TspI 


which is a combination of both, the incompressible Toner-Tu 


22[ 23 and the Swift-Hohenberg model 24 . The Swift-Hohenberg 


dynamic theory 
theory for flocks 

equation describes instabilities at finite wavelengths whereas the Toner-Tu theory is 
the basis for the description of flocking phenomena. A combination of both - a Tur¬ 
ing instability in the velocity field together with convective particle transport - is 
sufficient to explain the pattern formation observed in bacterial suspensions 121 . In 
particular, the experimentally observed velocity correlation function was reproduced 
by this hydrodynamic theory. However, the transport coefficients in the phenomeno¬ 
logical model remain unknown constants. 

In order to gain a deeper understanding of this pattern formation process, a par¬ 
ticle based model was suggested by us 25 . In this model, particles interact via com¬ 


peting velocity-alignment mechanisms: nearby particles tend to align their velocities 
whereas they anti-align their direction of motion with partic les which are further 
away. The idea of competing interactions goes back to Turing 26 who showed that 


instabilities on finite length scales arise in chemical systems due to the interplay of 
two competing species: a slowly diffusing activator and an inhibitor which diffuses fast 
compared to the activator. Within our active particle system, activation refers to the 
creation of local polar order by velocity-alignment which is inhibited by an opposing 
effect acting on larger length scales (anti-alignment). Similar types of competing inter¬ 
actions were studied in the context of spin models (ferromagnetic coupling to nearest 
neighbors vs. anti-ferromagnetic interaction with next-nearest neighbors 
the Kuramoto model for synchronization (conformists vs. contrarians 130 


27 


29 ) and 
However, 


note that - in contrast to these examples - the competing alignment interaction in our 
model leads to instabilities of the velocity field which is, in turn, related to convective 
mass transport resulting in numerous novel pattern formation phenomena 


25 


Generally, velocity-alignment interactions are understood as simple models for 
the interplay of various physical forces such as steric repulsion or hydrodynamics. 
The simplest active particle model with velocity-alignment interaction is the seminal 
Vicsek model [31 : particles move in continuous space at a constant speed and align 


their direction of motion with the local mean velocity plus some deviations (noise) 
in discrete time steps. Over the years, numerous modifications of the original Vicsek 
model were discussed (see e.g. [32]). In some studies, anti-alignment was included in 
different ways: for example, Menzel considered anti-alignment between two different 
particle populations 33 , whereas Lobaskin et al. studied a model where particles do 
only align if their velocity are already aligned to a certain degree 34 . For an overview 


on pattern formation in self-propelled particle systems, we refer the interested reader 


to a number of comprehensive review articles [2- 

In this work, we address the following questions. In section [2] we introduce our 
active particle model with competing alignment interaction in detail. In the subse¬ 
quent section, we report the results of numerical simulations and discuss the collec¬ 
tive dynamics exhibited by the system. Whereas we focused on the high density limit 
in [25] , we study the phase space in detail including low and high densities here. In 
section [4] we derive hydrodynamic equations in the mean-field limit allowing us to 
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understand the observed pattern formation qualitatively and to link our model to 
previously discussed phenomenological theories. Moreover, we estimate the quality of 
the approximations involved in the coarse-graining process and discuss their range of 
applicability. 


2 Self-propelled particle model 


We study an ensemble of N interacting self-propelled particles (SPPs) in two spa¬ 
tial dimensions. We assume that particles move at constant speed which is justified 
if all processes related to the self-propulsion mechanism occur on negligibly small 
timescales. The generic dynamics of interacting SPPs is determined by the following 
set of stochastic differential equation^ 


drj(t) 

d t 


f cos (f 7 
\sin <pi 


d <Pi(t) 
d t 


N 

^/2/W, (/) + ^ Pii Pj) 

i¥* 


(1) 


where we denote the spatial position of the iths particle at time t by rj(t) and its 
direction of motion by < Pi(t ). The first term on the right-hand side of the angular 
dynamics represents noise indicating the random reorientation of particles due to 
spatial inhomogeneities or fluctuations of the driving force, whereas the second term 
accounts for binary interactions among the particles. 

In the absence of inter-particle coupling, a SPP performs a persistent random 
walk: a particle moves ballistically at short timescales, whereas the motion is dif¬ 
fusive on large timescales due to angular fluctuations. The details of this diffusion 
process depend on the properties of the fluctuations 35 37 . Here, we consider the 
simplest case of Gaussian white noise with zero mean, (£,(£)) = 0, and temporal 6- 
correlations: = Sij 5(t — t'). The parameter D v parametrizes the intensity 

of the angular fluctuations. 

The binary interaction G(r ?; ;, ipj) generally depends on the relative coordinate 
r ji = Vj — i - ; of the particles i and j as well as on their individual orientations. Here, 


we consider the interaction scheme proposed recently by the authors in 25 including 
the following physical effects: 

(i) Due to the finite size of particles, these objects repel each other if they come 
closely together. 

=>• short-range repulsion 

(ii) In general, short-ranged anisotropic steric effects in combination with self-propulsion 
will lead to local orientational ordering. Here, we consider polar alignment: 

=>• velocity-alignment 

(iii) An active particle in a fluid or a polarly ordered cluster of swimming particles 
will induce hydrodynamic flows that in turn influence other particles. It has 
been shown for model swimmers that hydrodynamic interactions are generally 
anisotropic and may lead to velocity-alignment, anti-alignment or nematic align¬ 
ment depending on the relative orientation and position of two swimmers as well 
as the swimming type (see e.g. 38 ). Here, we replace the very complex hydro- 
dynamic interaction and take one additional effective interaction mechanism into 
account: 

=> anti-velocity-alignment. 


1 We work in natural units such that the speed, the mass of particles and the anti-alignment 
range are equal one. 
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Fig. 1. Visualization of geometric quantities and scheme of the pair-wise interaction. 


We parametrize the binary interaction in @>by 

U(r ji, <Pi, ipj) = «(|ryi|) sin (ipi - aji) + n(\rji\) sin (ipj - ipi). (2) 

The first term represents soft-core repulsion. The angle ctji = arg(rji) denotes the 
polar positional angle of particle j in the frame of reference of particle i (see also 
Fig. 0 for a graphical illustration). We assume a constant repulsion strength kq for 
|r_/i| < and vanishing repulsion otherwise: Ac(|rjj|) = Ko@(£ r — |rjj|). 

The second term describes the competition of alignment and anti-alignment: For 
/i(|rji |) > 0, particles tend to align their direction of motion whereas they anti-align 
in the opposite case. Throughout the paper, we will assume the following functional 
dependence of the alignment strength on distance 


Ml r ji|) = 


_ f d"M+ (l r ji|) 0 — l r ji| — £aj 

1~A 4 - (l r j*l) Cq ^ l r j'i| 5; 1) 


(3) 


where /i+ and are positive constants. The distance dependence of the competing 
interactions is determined by the functions A±(|rjj|). We normalize the distance de¬ 
pendence such that A± (|rj f |) £ [0,1]. Accordingly, fi + and denote the maximal 
alignment and anti-alignment strength, respectively (cf. Fig. |T|) . Following 25 
use a piecewise quadratic distance dependence in this study: 


we 


Wfrd-l-I^A A (r-)~ 4 ( r * ~ &X 1 ~ r *) 


(4) 


Note that the interaction described above is strictly short-ranged, but nonlocal. 


3 Simulation results - pattern formation and orientational order 

In this section, we report on simulations of the microscopic dynamics for various 
parameters which were performed to obtain a general understanding of the collective 
dynamics and pattern formation exhibited by the system. Our simulations revealed 














Will be inserted by the editor 


5 



PoM+ 


Fig. 2. The figure shows the part of the parameter space which we study in this work. 
Points indicate the parameter values for which simulation results are shown (cf. Fig. [d] - 
Fig. 10 1 . The background coloring (green, red, gray, white) was done in accordance with the 
predictions of the kinetic theory: in the green region, we expect the disordered (spatially 
homogeneous, isotropic) state to be stable; in the white area, we observe the emergence of 
polar order (local and global); in the gray and red colored part of the parameter space, 
we observe the emergence of patterns with a characteristic length scale (see section [d] for 
details). 


that the system exhibits various types of spatio-temporal dynamics which we will 
discuss in the following. 

We used the following fixed parameters in all simulations: repulsion range £ r = 0.1, 
repulsion strength k 0 = 10 and alignment range £ 0 = 0.2. We focus our discussion on 
the variation of the alignment and anti-alignment strengths. The parameter values 
which we are going to discuss are represented in Fig. [2] 

We go beyond our previous work 25 and explore not only the high density regime 
(po = 150) but also pattern formation at low densities (po = 20). The studied density 
values may misleadingly appear as extrem ely large in comparison to the values typi¬ 
cally studied in the Vicsek model 31 39 . However, the alignment range is less than 


one (£ q = 0.2) in our case. Hence, the mean number of particles within a box of edge 
length £ a is (n) = 0.8 for p 0 = 20 and (n) = 6 for p 0 = 150. This is comparable to 
the regimes typically studied in the context of the Vicsek-model. 

All simulations were performed using a rectangular domain with periodic bound¬ 
ary conditions with a linear size of L = 20 = 100 x £ a = 200 x £ r . This corresponds 
to a rather moderate system size with respect to the largest scale of interaction. We 
emphasize that the terms phase or state are therefore used here to describe stereo¬ 
typical patters exhibited by the system at the corresponding mesoscale and shall not 
be confused with phases in the thermodynamic limit. Accordingly, global order refers 
here to the ordering of the entire simulated system and does not imply long-range 
order. 


The initial conditions in all simulations were random positions and orientations 
of individual particles (disordered, homogeneous state). The simulation results were 
analyzed after the system had reached a stationary (dynamical) state, where no drift 
in order parameters and changes in patterns were observed anymore. The numerical 
time step was chosen sufficiently small to ensure numerical stability: At < 5 • 10~ 3 . 
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We measured polar and nematic orientational order using the corresponding order 
parameters 


*'=N 


N 

E 

3 = 1 






1 

N 


N 

Y e 2iv,j 
3=1 


(5) 


Furthermore, we calculated coarse-grained density p c (r,t) and velocity v c (r, t) by 
binning of the particles positions and velocities on a rectangular grid with the coarse- 
graining scale set by the linear dimension of a grid cell: dl = 0.2. Moreover, we 
calculated the momentum field via w c (r, t) = p c { r, f)v c (r, t). From the coarse grained 
fields, we deduced the two-point autocorrelation functions 


C pp (r) = p c (r',t)p c (r'+ r,t), C ww (r) = w c (r', t) • w c (r' + r, t) (6) 

as well as their corresponding two-dimensional Fourier transforms S p { k) = T[C PP \ 
and S w ( k) = J-[C WW \. In Eq. the overline indicates temporal averaging as well as 
averaging over the reference points r'. 


3.1 Low density case ( p 0 = 20) 


In general, no orientational or spatial order can emerge if noise is dominant in com¬ 
parison to the interaction strengths and, consequently, the system remains disordered, 
spatially homogeneous and isotropic. 

In the limit of vanishing anti-alignment (p- = 0), the system reduces to a pure 
alignment model and the large-scale dynamics corresponds to the Vicsek model with 
the usual order-disorder transition from a disordered to a polarly ordered state at 
high alignment and small noise. In particular, we observe the emergence of large-scale 
travellin g p olar bands close to the order-disorder transition typical for Vicsek-type 
models 


These structures can be clearly identified by the asymmetry of C t 


40 


pp 


and C ww reflecting an increased density and velocity autocorrelation perpendicular 
to the polar orientation and low-wavenumber peaks along the polar direction in both 
spectra S p>w (Fig. 41). Far away from the order-disorder transition (large p+ and low 
p~), the model exhibits the spa tial ly homogeneous phase with polar order described 
and analyzed by Toner and Tu 


22 


Small levels of anti-alignment do not destroy polar order at the system sizes stud¬ 
ied. Only if p- is increased beyond a critical value, we observe the breakdown of global 
polar order. At low density, this breakdown is accompanied by the formation of small 
polar clusters, which are unstable and frequently break up. The corresponding density 
autocorrelation and its Fourier transform are rotationally symmetric and decay fast 
from a central peak indicating no long-ranged spatial structure. However, the ringlike 
structure of C ww and the maximum of S w (k) at a finite |k| identify a typical length 
scale of velocity correlations corresponding to the typical size of clusters (Fig. pi. 

If the dynamics is dominated by anti-alignment (large p -), we observe the forma¬ 
tion of patches with nematic order where the polar clusters form lanes with neigh¬ 
boring lanes moving into opposite directions leading to global nematic order. This 
anti-parallel movement distinguishes these patters from the active crystals described 


41 


42 


C pp and S p show a high degree of spatial structure corresponding to the 
characteristic distance between clusters and lanes, respectively. Interestingly, the mo¬ 
mentum autocorrelation function C ww also shows crystal like structure corresponding 
to a hexagonal lattice (Fig. pi). We emphasize that the autocorrelation functions were 
obtained from time-averagea density and momentum fields. In general, these autocor¬ 
relation functions and spectra will differ from ones obtained from a single snapshot 
at a fixed time t. 
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Surprisingly, the state with global nematic order and local polar order also emerges 
for anti-alignment only (/.t + = 0) as shown in Fig. [b] In this case, particles do not 
interact for distances £ r < rji < £ a and local polar order is enforced by the strong 
anti-alignment at larger distance. The formation of the particular dynamic cluster 
state with dense polar clusters moving in lanes is the only possibility to ensure a 
stationary pattern when anti-alignment interaction is dominant. Nematic lanes were 
reported previously from other active matter systems (see e.g. 43-45 and |[5] for a 


recent review as well as references therein), however not in absence "of an (effective) 
alignment interaction and never with additional periodic order within individual lanes. 


3.2 High density case ( p 0 = 150) 


We discussed simulation results in the high-density regime in 25 . Here, we review 


these previous observations and go beyond by exploring a wider range of parameters. 
Snapshots of simulations in the high density regime are shown in Fig. [3] 

Similar to the low-density case, we observe the emergence of large scale polar 
order for vanishing /r_. However, there are no bands due to the high density and we 
observe the Toner-Tu phase (Fig. [7]). An increase in /i_ beyond a critical value leads 
to the breakdown of polar order. Just above this transition, we observe a dynamical 
regime without any spatial order but with a finite scale in the momentum field, 
which is reflected by a single ring of negative correlation about the central maximum 
of C ww and the rotationally symmetric S w with a maximum at finite |k|. In contrast 
to the low-density case, we observe extended convective flows (jets) and transient 
vortices (Fig. [8]) instead of diffuse clusters. This is the regime discussed in 25 which 
reproduces qualitatively and quantitatively the observations of mesoscale turbulence 
in dense bacterial suspensions. 

An increase in initially leads to the stabilization of the vortices and we observe 
the emergence of patches of spatially ordered vortex arrays arranged in a rectangular 
lattice with rather homogeneous density (Fig. [9]). However, a further increase in /i_ 
leads also to growth of density inhomogeneities and eventually the rectangular vortex 
pattern transforms into a hexagonal array of rotating clusters (Fig. 10 ). 



Fig. 3. Close-up snapshots in the high density regime. From left to right: homogeneous 
polarly ordered state, mesoscale turbulence, rectangular vortex array, rotating clusters. 
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Fig. 4. Polar bands ($ v = 0.65 and = 0.33) at strong alignment p+ = 4, weak anti¬ 
alignment p- = 0.004, low density po = 20 and noise D v = 1.0. (A) Global snapshot with 
particle positions shown by symbols and orientation indicated by the color according to the 
inset. (B) Zoom into a subregion of size 5x5 with particle positions and orientations shown 
as vectors. (C) Time averaged density autocorrelation function C pp ( r) and (D) its Fourier 
transform S p ( k). (E) Time averaged momentum field autocorrelation function C ww ( r) and 
(F) its Fourier transform S w ( k). 



Fig. 5. Disordered clusters (d> p = 0.021 and = 0.020) at strong alignment p+ = 4, 
moderate anti-alignment p~ =0.1, low density po = 20 and noise D v = 1.0. (A) Global 
snapshot with particle positions shown by symbols and orientation indicated by the color 
according to the inset. (B) Zoom into a subregion of size 5x5 with particle positions and 
orientations shown as vectors. (C) Time averaged density autocorrelation function C pp ( r) 
and (D) its Fourier transform S p ( k). (E) Time averaged momentum field autocorrelation 
function C ww ( r) and (F) its Fourier transform S v ( k). 
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Fig. 6. Nematic order (<P p = 0 and d? n = 0.93) due to lanes of polar clusters without 
alignment p+ = 0 at strong anti-alignment p_ = 4, density po = 20 and noise D v = 1.0. 
(A) Global snapshot with particle positions shown by symbols and orientation indicated 
by the color according to the inset. (B) Zoom into a subregion of size 5x5 with particle 
positions and orientations shown as vectors. (C) Time averaged density autocorrelation 
function C pp ( r) and (D) its Fourier transform S p ( k). (E) Time averaged momentum field 
autocorrelation function C ww ( r) and (F) its Fourier transform S„(k). 



Fig. 7. Large-scale polar order (4> p = 0.92 and <L> n = 0.73) at high density po = 150 
with p+ = 1.6, p_ = 0.0016 and D v = 1.0. (A) Global snapshot of the density field. (B) 
Global snapshot of the orientation field with polar orientation indicated by color according 
to the inset. (C) Time averaged density autocorrelation function C pp ( r) and (D) its Fourier 
transform S p ( k). (E) Time averaged momentum field autocorrelation function C ww ( r) and 
(F) its Fourier transform S w ( k). 
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Fig. 8. Mesoscale turbulence - complex spatio-temporal convective flows = 0.01 and 
<&„ = 0.01) - at high density po = 150 with p+ = 1.6, p~ = 0.008 and D v = 1.0 (see 
also [25] ). (A) Global snapshot of the density field. (B) Global snapshot of the orientation 
field with polar orientation indicated by color according to the inset. The black-white elipse 
indicates a jet and the ring a transient vortex. (C) Time averaged density autocorrelation 
function C pp ( r) and (D) its Fourier transform S p ( k). (E) Time averaged momentum field 
autocorrelation function C ww ( r) and (F) its Fourier transform S w ( k). 



Fig. 9. Vortex lattice = 0.0 and = 0.01) with (local) rectangular positional order at 
high density po = 150 with p+ = 3.2, p- = 0.032 and D v = 1.0 (see also |25| ). (A) Global 
snapshot of the density field. (B) Global snapshot of the orientation field with polar orien¬ 
tation indicated by color according to the inset. (C) Time averaged density autocorrelation 
function C pp ( r) and (D) its Fourier transform S p ( k). (E) Time averaged momentum field 
autocorrelation function Cu^r) and (F) its Fourier transform S ro (k). 
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Fig. 10. Dense rotating clusters = 0.06 and <?„ = 0.02) at high density po = 150 
with p,+ = 1.6, g~ = 1.024 and D p = 1.0. The velocity field as well as the density are 
periodically modulated. (A) Global snapshot of the density field. (B) Global snapshot of 
the orientation field with polar orientation indicated by color according to the inset. (C) 
Time averaged density autocorrelation function C pp ( r) and (D) its Fourier transform S p ( k). 
(E) Time averaged momentum field autocorrelation function C ww ( r) and (F) its Fourier 
transform S w ( k). 


4 Kinetic & hydrodynamic theory 

In this section, we turn our attention to the analytical analysis of the microscopic 
model aiming at a qualitative understanding of the mechanism leading to pattern 
formation on macroscopic scales. For that purpose, we systematically derive kinetic 
as well as hydrodynamic equations which describe the dynamics of the system in this 
regime. In particular, we quantify the quality of the approximations involved in the 
coarse-graining process and discuss their range of applicability. Our analysis is based 
on the study of the one-particle distribution function 

N 

= ^^ (2) (r-r l (f))(5(v?-^ i (t))^), (7) 

i—1 

which determines the density of particles at position r moving into the direction 
<p. Its dynamics cannot be obtained rigorously in a closed form for an interacting 
many-particle system. This fact is easily understood: In order to calculate the force 
between two particles, one has to know the probability to find a particle at position 
r given that another particle is located at r'. Hence, the dynamics of the one-particle 
distribution involves the two-particle probability distribution. Here, we overcome this 
hierarchy problem by assuming that fluctuations are small ( mean-field assumption) 
enabling us to approximate the two-particle probability distribution by the product 
of one-particle densities (see appendix [A] for details). Hence, our kinetic theory does 
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not account for inter-particle correlations. However, this assumption is justified in the 
limit of high particle densities: In this situation, relative density fluctuations are small 
and, consequently, every particle interacts with an approximately constant number 
of neighbors 


25 46-48 


Using the mean-field arguments above, we obtain the following nonlinear Fokker- 
Planck equation for the dynamics of the one-particle density 


shr p) ' Vr ^ r ’ ip { F[p] P(r ’ V’ ^ 


dp(r,ip,t) 

dt 


where the nonlinearity enters through the force 


dp( r, ip, t) 
dip 


( 8 ) 


F\p] = jd 2 r'J d^'{K(|r'|)sin[^-arg(r')]+ir(|r'|)sin[^' - ip]}p(r+r',ip',t). (9) 


The terms in Eq. ([8]) arise from the following microscopic rules: The first term de¬ 
scribes convection due to active motion at constant speed, the second term reflects 
the binary, nonlocal interaction (in mean-field approximation) and the third term 
accounts for angular fluctuations. This nonlinear Fokker-Planck equation is solved by 
the spatially homogeneous, rotationally invariant state 


Po = 


Po 

27T 


( 10 ) 


which describes a disordered configuration without orientational order and a spatially 
homogeneous particle density po- In the following paragraphs, we study the linear sta¬ 
bility of this fixed point considering the full dynamics (I8| and compare the results to 
the predictions of the linear stability analysis using a reduced set of equations ( hydro- 
dynamic theory). This comparison allows a critical assessment of the approximations 
used in the derivation of the hydrodynamic theory. 


4.1 One-particle Fokker-Planck equation in Fourier space 

In order to study the nonlinear Fokker-Planck equation ([8]) analytically, it is conve¬ 
nient to work in Fourier space with respect to the angular variable <p. Here, we use 
the following convention for the Fourier transformation: 


fn{ r,t) 


r 2n 


dip e mv p(r, ip, t). 


( 11 ) 


Due to the nonlocality of the Fokker-Planck equation - the dynamics of p(r,ip,t) 
depends on the integral over p(r, < p,t) itself - the Fourier transformation is somewhat 
involved. It is helpful to express the force F[p], Eq. Q, in terms of Fourier modes 
first. For this purpose, we assume that the one-particle density p( r + r', ip,t) can be 
expanded in a Taylor seriesj^] in the spatial coordinate r' around r' = 0. In other 
words, we transform the integral operator in a series of differential operators. We 
keep all orders in the Taylor expansion and obtain 


F = 


(/Wi(r,t) - /CAV/ 0 (r,f)) — c.c 


(12) 


2 


Furthermore, we assume the global convergence of the Taylor series. 
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tors 


In (121, we introduced the complex derivative V = ^ ( d x + id v ) as well as the opera- 


HA = 


n =0 


and 


fc-A = 2 / dr re(r) 


dr r/i(r)Jo {ry/—Aj = 


drr p{r) 




= £ 

71 =0 


dr r ft(r) 


4”(?r!) 2 


2n+l 




4"n!(n + 1)! 


A* 


(13) 


(14) 


which stem from the Taylor expansion. We denote the Laplace operator by A. The 
first term in (131 and (141 is written in a formal way using Bessel functions of the first 
kind J„{x), whereas the second follows immediately from the series expansion of the 
Bessel functions. These series reflect the nonlocality of the interaction. Interestingly, 
the alignment interaction couples to the first Fourier mode only (proportional to 
the polar order parameter) and the repulsive interaction couples to the zeroth mode 
(particle density). 


Having rewritten the force F\p] as it is done in Eq. (121, the dynamics of the 
angular Fourier coefficients is directly obtained by taking the temporal derivative of 
(111 and inserting the nonlinear Fokker-Planck equation for the one-particle density 
on the right hand side: 


dfn 

dt 


= - (y/„_i + V7„ +1 ) + nir 


fn -1 fi-A f 1 — fn+1 &A f -1 


fn— 1 KjV/o - fn+1 JC A V*fo - n 2 D v f n . 


(15) 


The dynamics of the nth Fourier coefficient (151 is in general coupled to others 
meaning that the nonlinear, nonlocal Fokker-Planck equation is equivalent to an in¬ 
finite hierarchy of equations in Fourier space. This hierarchy is, however, insightful 
because it allows the investigation of the consequences of approximations, such as the 

— negligence of Fourier coefficients f n for \n\ > n c , where n c > 0, 

— negligence of derivatives of a particular order, 

which are studied within linear stability analysis of the spatially homogeneous state 
in the following paragraphs. 


4.2 Linear stability analysis within kinetic theory 

Before we discuss the consequen ces of several approximations, we illustrate the pre¬ 
dictions of the full kinetic theory (15) in this section. In particular, we study the linear 


stability of the disordered state by linearizing the kinetic equations for small pertur¬ 
bations. In Fourier space, the dynamics of the perturbations Sf n fk,t) is determined 
by a linear system of coupled differential equations 




(16) 


l= — C 


with the coefficients 


M ni ( |k|) =- [(k x + ik v )S n - 1> i + (k x - ik y )5 n+u } + mrp 0 fi k (S nA - 5 n _i) 5„,z 

^ D v 5 n ,i -\- ~ /Cfc [(&£ T ik y )5 n + {kx iky)5 n .—i] 5p o • (17) 


2 
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In ( [l7| ), the Fourier representation (spectra) of the nonlocal operators fiA and Ka 
were introduced: 

Afc = J drrn(r)J 0 (\k\r), JC k = 2 J dr rn{r) ■ ( 18 ) 


The eigenvalues of the matrix M as a function of the wavevector k - also referred 
to as dispersion relation - determine the stability of the disordered state. If the real 
part of an eigenvalue is positive, the disordered state is unstable. Since the disordered 
state is invariant under spatial rotations, the dispersion relation does depend on the 
wavenumber |k| only. We denote the dispersion relation by cr(|k|) and its real part by 
5J[er(|k|)]. Once cr(|k|) has been calculated, we can distinguish different scenarios: 

(i) If 3i[<T(|k|)] is positive at its maximum located at |k| = 0, the spatially homoge¬ 
neous system is unstable. 

(ii) If R[<r(|k|)] is positive at its maximum located at a finite value of |k|, a pattern 
with a finite characteristic length scale is expected to emerge (Turing instability). 

(iii) If 5R[cr(|k|)] is non-positive, the disordered state is stable. 

Information about case (i) is obtained from the study of the spatially homogeneous 
system by setting |k| =0. The matrix © reduces to a diagonal matrix in this case. 
Hence, the eigenvalues are trivially obtained. The analysis reveals that only the first 
Fourier coefficient (related to the polar order parameter) can become unstable. This 
instability is due to the alignment interaction which does only couple to the first 
Fourier coefficient as can be seen from Eq. ( fl2| ). The parameter region where this 
instability occurs - signaling the emergence of local polar order - is determined by 


npo 


drr/j(r) ) - D v > 0. 


(19) 


Equation (19I is a generalization of the expression derived in 46 for an arbitrary 


distance dependent velocity-alignment strength. Consistently, spatially homogeneous, 
polarly ordered states cannot emerge if the interaction is dominated by anti-alignment, 
i.e. if the condition f 0 °° drr/j(r) < 0 is fulfilled. 

The full dispersion relation cannot be obtained analytically. In order to study 
the linear stability of the disordered state with respect to spatially inhomogeneous 
perturbations, we calculated the dispersion relation numerically. Since the matrix 
M is infinite dimensional, it must be truncated first. However, we can exploit that 
the Fourier modes /„ are damped proportional to n 2 due to angular fluctuations. 
Consequently, high order Fourier modes are strongly damped. Accordingly, we observe 
that the largest eigenvalues of the matrix © do not depend on the truncation used 
if a reasonably high number of Fourier modes is considered. Here, we reduce M to a 


(101 x 101)-matrix. The results of the linear stability analysis are shown in Fig. 11 


The linear stability analysis reveals the following regimes: If the alignment-strength 
and the density are small, the disordered state is linearly stable (parameter region 
(A) in Fig. 11). An instability at finite wavelengths (Turing instability) appears if the 
anti-alignment strength is high (domain (B)). This instability signals the transition 
towards an ordered state which possesses a characteristic length scale (cf. numerical 
simulation in section [3]) . The disordered state is unstable with respect to spatially 
homogeneous perturbations if the alignment strength is high and anti-alignment is 
negligible (domain (D)). In between region (B) and (D), an intermediate regime (C) 
exists where the dispersion relation is positive at |k| = 0 but possesses a maximum 
at a nonzero wavenumber. 

The occurrence of a Turing instability in our model is due to the concurrence 
of short-range alignment and anti-alignment interaction. This instability can only 
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PoM+ 



Fig. 11. Results of the numerical linear stability analysis of the disordered state po, Eq. (110]) , 
within kinetic theory. Top: The (pop+)-(p-/p+)-plane of the parameter space. Bottom: 
Exemplary dispersion relations K[cr(|k|)] obtained numerically from (17l. If the dispersion 
relation is positive, the disordered state undergoes an instability. In parameter region (A), 
the disordered state is stable. We find Turing instabilities in the parameter range ( B ). 
In (C), the dispersion relation possesses a maximum at finite |k| and, additionally, the 
homogeneous system (|k| = 0) is unstable. In region (D), the spatially homogeneous system 
is unstable. The boundary cr( |k| = 0) = 0, determined by (191, is represented by a dashed 
line. Arrows in top panel indicate the shift of critical lines when the noise intensity D v is 
increased. Parameters of dispersion relations: poP+ = 1> P-/p+ = 10” 1 in (A); poP+ = 10 2 , 
p-/p + = 5- 10~ 2 in (B); pop+ = 10 2 , p-/p + = 10~ 2 in (C); pap+ = 10 2 , p-/p + = lCF 2 in 
(D); = 1., £ a = 0.2, = 0.1, «: 0 po = 1-5 ■ 10 3 . 
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be found if the anti-alignment interaction is sufficiently strong. Consequently, this 
instability mechanism is not present in models with purely aligning interactions. 

From the linear stability analysis, we cannot rigorously infer the structure of 
the emergent pattern which crucially depends on nonlinear terms. However, we may 
expect the emergence of collective states which possess a characteristic length scale 
in region (B) and (C). Indeed, we observe these state such as clusters of characteristic 
sizes, mesoscale turbulence and vortex arrays in Langevin simulations of the particle 
dynamics (cf. section [3] in particular Fig. |2j) . 


4.3 Hydrodynamic theory 


The linear stability analysis of the disordered state within kinetic theory revealed the 
existence of several instability mechanism. However, the dispersion relations could 
only be obtained numerically. In this section, we discuss the possibility to describe 
the instabilities of the disordered state analytically by a reduced set of linearized 
hydrodynamic equations aiming at a qualitative understanding of the macroscopic 
pattern formation process. 

The basic assumption of the hydrodynamic theory is that only those observables 
in the hierarchy (15) are important whose large scale dynamics is slow either due to 
the fact that they describe a conserved quantity or they describe broken continuous 
symmetries. For the identification of the relevant degrees of freedom in our model, 
we exploit that (i) the particle number is conserved and (ii) the microscopic inter¬ 
action favors local polar order due to velocity-alignment. Hence, the relevant order 
parameters are the particle density 


P(r,f) = /o(r,t) = / d<pp(r,<p,f) 
Jo 


( 20 ) 


and the momentum field 


w (r , f ) = l d^(r , V ,t) = (KJ) (21a) 

which is proportional to the local flow field v(r,t) = w(r,f)/p(r,f). Please note that 
we consider the order parameters locally in space. The Langevin simulations of section 
[3] revealed that all macroscopic patterns display local polar order. This assumption 
does not exclude the description of nematic states such as polar clusters moving in 
opposite directions (cf. Fig. [b]) which possess local polar and global nematic order. 

In general, the dynamics of the density and the momentum field is coupled to 
higher order Fourier modes. In particular, the dynamics of the polar and nematic 
order parameter is not independent. However, we argued above that the dynamics 
of the relevant order parameters is slow on large scales. The assumption of fast re¬ 
laxation of higher order Fourier modes enables us to truncate the hierarchy (15) at 


low order. In particular, we express the nematic order parameter as a function of 
the polar order parameter and the density by adiabatic elimination, and further ne¬ 
glect all higher order Fourier coefficients following the same reasoning. Here, we focus 
on linear terms since we are interested in comparing the predictions of the linear 
stability analysis of the disordered state. We obtain the following set of linearized 
hydrodynamic equations: 


— = — V • w, 


dp 

dt 

dw 

~dt 


1 npoKog 

2 + 6 


Vp - D v w + npopAW + 


Aw 
16 A. 


(22a) 

(22b) 
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Fig. 12. Comparison of the results of the linear stability analysis of the spatially homoge¬ 
neous state using the Fokker-Planck equation |8} and the hydrodynamic approximation (22 I. 
The results corresponding to the kinetic theory are shown in the same way as in Fig. Ill 
Predictions of the hydrodynamic theory are shown by lines. In the high noise limit (right), 
the hydrodynamic approximation is exact, whereas deviations are found for low noise values 
(left). Parameters: D v = 0.1 (left), D v = 10 (right). Other parameters: see Fig. 
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The continuity equation (22a I reflects the conservation of the number of particles. 
The equation for the momentum field contains several physical effects. The first term 
which is proportional to the gradient of the density describes a macroscopic particle 
flow compensating density gradients. The prefactor in brackets is proportional to the 
pressure which contains two contributions: (i) the constant term is due to the active 
motion; (ii) the second term is proportional to the repulsive force which effect ively 
increases the pressure, i.e. it reduces the compressibility. The second term in (22b I 
accounts for the effect of noise counteracting the emergence of polar order. The non¬ 
local alignment inter act ion is described by the operator p, z \ defined via a series of 
Laplacian operators (13). The last term represents the reduction of the effective vis¬ 
cosity due to angular noise. It is the only term which is influenced by the truncation 
scheme used to derive ([22 1. 


Linear stability analysis within the hydrodynamic theory 


The truncation of the infinite hierarchy of Fourier modes to 
hydrodynamic equations is a considerable simplification. In 
analytical calculation of the dispersion relations which can 
Fourier transformation of the hydrodynamic equations (221 
diagonalization: 


the simplified system of 
particular, it allows the 
be directly obtained by 
and subsequent matrix 


C7i(|k|) = npoflk - D v 


02,3(|k|) 


°i(|k|) ± 


lk| 2 

16 iV 

]J M|k|)] 2 - 2 ^1 + Tr ^ptl^J |k| 2 


(23a) 

(23b) 


We address the question in which parameter range this approximation is applicable 
and whether these approximated dispersion relations can be used to describe the 
destabilization of the disordered state. 

The results of the parameter study are illustrated in Fig. [12] The analysis was 
performed for high and low noise values D For high noise amplitudes, we obtain 
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perfect agreement of the predictions of the hydrodynamic equations (22I with the 
results of the linear stability analysis which was based on the one-particle Fokker- 
Planck equation ([8]). However, the predictions of the hydrodynamic theory deviate 
from the kinetic theory in the limit of low noise values. These deviations are due to 
the adiabatic elimination of the second Fourier mode and the negligence of all higher 
Fourier coefficients. These deviations could have been expected: The assumption that 
high order Fourier modes relax fast is only justifiable for high noise values, since D^ 
determines the relaxation rate of the Fourier coefficients as can be seen from (15I . 

Consequently, we conclude that it is possible to describe the linear dynamics of the 
system by the hydrodynamic equations (221 as long as the noise strength is high, i.e. 
a considerable timescale separation between the lowest and higher Fourier coefficients 
exists. 


Approximating the nonlocal alignment interaction 


As already mentioned in section |4.1| one can possibly approximate the infinite series 
fiA of Laplacian operators by a truncated polynomial: 


N 0 


V-A ~ PA,N C 


= E 

n—0 


dr rp(r) 


4 ™(n!) 2 


A n . 


(24) 


Naturally, this raises the question how to choose N c (order of truncation) and how 
reliable the corresponding approximations are. The stability of the dynamics requires 
that 


(- 1 ) 


JV C 


dr p,(r) 


„2N C +1 


Po 


< 0 


(25) 


holds. As a rule of thumb, (25) reflects that N c must be odd if alignment interaction 


dominates and N c must be even if anti-alignment is predominant. 

What is the nature of an approximation like p41)? The operator p,A describing 
the nonlocal (anti-)alignment interaction is linear and, thus, it possesses a character¬ 
istic spectrum introduced in Eq. (18). Approximating the operator by a series of 
Laplacian operators, cf. ( [24] ), is equivalent to approximate the spectrum by a Taylor 
series in powers of the wavenumber |k| around |k| = 0 restricting the applicability of 
this approximation to small wavenumbers (large length scales). The stability require¬ 
ment (251 ensures that the spectrum at large wavenumbers is negative. Otherwise, 


fluctuations on small length-scales would be amplified without any bound. 

On phenomenological grounds, there are two conceivable types of spectra: 

(i) the spectrum possesses a maximum at |k| = 0 and decreases monotonically; 

(ii) the spectrum possesses a maximum at a finite wavenumber. 

The former case was studied in the context of the Vicsek model (pure alignment 
interaction 31 ) in the corresponding field theory by Toner and Tu 49]. The latter 
was recently proposed by Dunkel et al. 50 in order to describe mesoscale turbulence. 


Using the expansion of the operator pAi we can map our linearized hydrodynamic 
theory to these phenomenological equations 


dp 

dt 

dw 

~dt 


— V • w, 

— ^ + 7rp ° l ^°^ r ^ Vp + aw + r 0 Z\w + r 2 A 2 w, 


(26a) 

(26b) 
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Fig. 13. Domains of stability and instability of the disordered state. Predictions of the 
kinetic theory ([8} are depicted as gray shaded regions (cf. Fig. 111. The colored lines are 
derived from the dispersion relations of the linearized hydrodynamic theory (22 l for different 
values of the truncation parameter N c , cf. Eq. (24 1 . The critical line of the Turing instability 
(blue dashed line separating (A) and (B)) is not correctly predicted in the case N c = 2. 
Qualitative agreement is achieved for N c = 4, higher order approximations are almost exact. 
Parameters: D v = 10, others as in Fig. 


and derive the dependency of transport coefficients on microscopic model parameters 


r 0 

r 2 


TTPo 

TTPo 

4 

Trio 

64 



(27a) 

(27b) 

(27c) 


50 , i.e. only if the anti-alignment 


The linearized hydrodynamic equations can describe a Turing instability of the disor¬ 
dered state if the effective viscosity r 0 is negative 
interaction is sufficiently strong. 

It is a priori not clear how many terms are necessary in order to reproduce quanti¬ 
tatively the instabilities found from kinetic theory using a truncated operator £ia,n c - 
We close this paragraph by answering this question for our model. In the case of 
N c = 0, the results deduced from the spatially homogeneous equations are recovered 
(cf. Eq. (19 >). The results of the linear stability analysis for N c G {2,4, 6, 8} are shown 
in Fig. m It is evident from our analysis that it is insufficient to approximate (ia in 
the leading order N c = 2. The disagreement is due to the fact that not all relevant 
length scales are resolved properly in this case. However, the situation is essentially 
different if the approximation N c = 4 is used. This approximation enables the qual¬ 
itative prediction of the structure of the phase space. Higher order approximations 
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N c > 6 can be considered as almost exact on the linear level. We stress that the num¬ 
ber of terms which are necessary to reproduce the spectrum sufficiently well, i.e. 
how many terms are necessary to resolve all relevant length scales, is not a universal 
number but depends on model details. 


4.4 Conclusions 


The analysis presented in this chapter reveals the existence of a Turing instability of 
the coarse-grained velocity field due to the anti-alignment interaction - there exists 
a region in parameter space where the dynamics of the system is linearly unstable at 
finite wavelengths. Furthermore, we showed that a system of hydrodynamic equations 
(261, where the nonlocal operator jx^ is approximated by a truncated series, provides 
an excellent approximation for the prediction of phase boundaries given the following 
conditions are fulfilled: 


— The dynamics of the Fourier coefficients decouple, i.e. the truncation of the hierar¬ 
chy of Fourier modes is possible, if a considerable timescale separation among the 
modes exists. In particular, high Fourier modes are strongly damped if angular 
fluctuations are large justifying the negligence of these high order modes. This 
property allows to truncate the infinite hierarchy of Fourier modes and enables us 
to replace it by a low dimensional system. 

— The nonlocal interaction operator must be approximated sufficiently well, i.e. the 
series expansion of fiA must contain enough terms in order to resolve all relevant 
length scales. In our particular model, this is the case if all terms including the 
term proportional to Z\ 6 w are taken into account. 


The linearized theory discussed here is not capable of predicting the macroscopic 
patterns since nonlinear terms were not considered for simplicity. However, the the¬ 
ory confirms the existence of a Turing instability of the polar order parameter in 
our model. Such an instability mechanism was proposed in 50 on phenomenological 


grounds and sh own to be sufficient to explain mesoscale turbulence in dense bacterial 
In this work, we provide a microscopic explanation of this insta- 
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suspensions 

bility: the concurrence of short-range creation of local polar order with an opposing 
effect on larger scales. 


5 Summary 


We studied the pattern formation in an ensemble of self-propelled particles which 
move at constant speed and interact via competing velocity-alignment interactions: 
at short length scales, particles tend to align their velocities whereas they anti-align 
their direction of motion with particles which are further away 25 . The interplay 


of alignment and anti-alignment interaction has to be understood as an effective 
description of two competing mechanisms: locally, polar order arises due to anisotropic 
steric interactions or inelastic collisions for example, whereas the emergence of global 
polar order is suppressed by an opposing effect, e.g. hydrodynamic interaction. 

Numerical simulations of the particle dynamics reveal the emergence of surpris¬ 
ingly complex spatio-temporal pattern formation. On the one hand, our model reduces 
to a continuous time version of the Vicsek model 131,46 for vanishing anti-alignment 


39 


interaction. Therefore, we observe a phenomenology similar to the Vicsek model 
in this limiting case. Especially, large scale traveling bands emerge in the low den¬ 
sity regime close to the onset of polar order. Moreover, numerous additional phases 
were found due to the anti-alignment interaction. If this interaction is strong and 
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the particle density is high, we observe vortex arrays and synchronously rotating 
cluster. By gradually lowering the anti-alignment interaction, the velocity-alignment 
acting on short length scales induces local polar order and, thus, mesoscopic con¬ 
vective flows. In this intermediate regime, the anti-alignment interaction inhibits the 
emergence of global polar order and, consequently, we observe a rotationally invariant 
state possessing local polar order but global disorder (mesoscale turbulence). Another 
characteristic of this state is the emergence of transient vortices, which are separated 
by a characteristic length. Moreover, we observe extended convective flows (jets). At 
low densities, the actual patterns differ strongly: We observe small polar clusters mov¬ 
ing randomly for intermediate anti-alignment strength. At strong anti-alignment, the 
system exhibits large scale nematic order as polar clusters start to form lanes moving 
in opposite directions. 

Common to all the patterns is a characteristic length scale which manifests itself 
mainly in the momentum field. Studying the linear stability of the disordered state 
within mean-field theory, we were able to relate this pattern formation process to a 
Turing instability. Furthermore, we analyzed in detail the possibility of describing this 
instability using a reduced set of hydrodynamic equations. We showed that a system 
of hydrodynamic equations can quantitatively reproduce the relevant instabilities if 
two conditions are fulfilled: (i) the reduction to a low dimensional system of order 
parameter equations is only possible if a considerable timescale separation between 
fast and slow modes exists which is particularly fulfilled for high noise in our model; 
(ii) the nonlocal interaction has to be described sufficiently well by taking into account 
a high number of derivatives in the coarse-grained description. 

Based on the analysis of our self-propelled particle model, we suggest that the 
interplay of competing alignment interactions may generally give rise to a Turing 
instability in the polar order parameter field and that this mechanism in combination 
with convective particle transport can explain the emergence of mesoscale turbulence 
and vortex arrays in active matter systems. 
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A Mean-field approximation 

The derivation of the nonlinear Fokker-Planck equation ([8]) relies on a mean-field 
approximation (MFA) - the negligence of fluctuations - which we illustrate in detail 
in this section. The central quantity is the microscopic particle density V(r,ip,t) 
defined by 


N 


= si2) ( r ~ “ <?*(*)) 


(28) 


which is a stochastic fluctuating field variable. From this field, we can compute the 
one-particle probability density function to find a particle at position r with orienta¬ 
tion ip at time t via 
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as well as the two-particle probability density function to find one particle at position 
r with orientation p and another particle at position r 1 with orientation p' at time t 
via 


P 2 (r,<p-,r',tp'-,t) 


(E^ (2) ( r - r i (t))(5( 2 )(r' - Yj(t))5(ip - Pi{i))5{p’ - <Pj(t)j} 

N(N - 1) 

(30a) 

(V(r,ip,t)V(r\ip',t)) _ P 1 (r,p,t)5^ 2 \r - r')S(p - p') 

N(N-l) N-l ’ 


In mean-field approximation, fluctuating variables are approximated by their average 
values thus neglecting fluctuations: 


V(r,p,t) = (P(r, p,t)) + SV( r , p, t ) « (V(r, tp,t)). 

(MFA) 


(31) 


Therefore, we deduce from (30b) that we can approximate the two-particle probability 
density function by the product of one-particle probability densities: 

NP 1 (r,p,t)P 1 (r',p',t) Pi(r, p, t)5^ (r - r')5(p - p') 


P 2 (r,p; r',p ; t) 


(MFA) 


N-l 


N-l 


(32a) 

w P 1 (r,p,t)P 1 (r',p',t). (32b) 

(AT->co) 

Thus, the mean-field approximation is equivalent to the negligence of inter-particle 
correlations. 
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